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COSMIC RAYS IN A RANDOM MAGNETIC FIELD: 
BREAKDOWN OF THE QUASILINEAR DERIVATION 
OF THE KINETIC EQUATION 
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Abstract 

We consider the problem of deriving a kinetic 
equation for the cosmic ray distribution function in a 
rarriom magnetic field. A model is adopted which is 
mathematically simple but which contains the essential 
physics. We investigate the perturbation expansion upon 
which the quasilinear treatment employed by previous 
authors is based. As pointed out by Klimas and Sandri, 
the existence of resonant particles causes the breakdown 
of the adiabatic approximation frequently used in this 
theory. We find further that resonant particles cause 
a general secular growth of higher order terms in the 
expansion which invalidates the entire perturbative approach. 


*This work was supported by the National Aeronautics and 
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COSMIC RAYS IN A RANDOM MAGNETIC FIELD: BREAKDOWN OF 


THE QUASILINEAR DERIVATION OF THE KINETIC EQUATION 
I. Introduction. 

During the past few years several authors have considered 
the problem of how best to describe the evolution of a 
distribution of cosmic ray particles in a random . electromagnetic 
field. Certain of these investigations have employed a 
quasilinear analysis (Vedenov, Velikhov and Sagdeev 1961, 

1962; Drummond and Pines 1962) of the Vlasov equation to derive 
a kinetic equation of the Fokker-Planck form (Hall and Sturrock 
1967; Hasselmann and Wibberenz 1968; Kulsrud and Pearce 
1969; Jokipii 1971, 1972). This approach treats the effect of 
the random field as a perturbation of the orbits of particles 
moving in an average background field. The distribution 
function is expressed as a power series in a small parameter 
which characterizes the strength of the random field. 

Quasilinear theory assumes that terms in this series of higher 
than second order are negligible. In order to obtain a 
Fokker-Planck equation it is further assumed that there are 
two time scales: the correlation time of the fluctuating 
field as seen by a particle moving along its unperturbed 
trajectory, and the much longer relaxation time of the 
distribution function. 

Recently, these treatments have been criticized by 
Klimas and S&ndri (1971) on the grounds that there are cosmic 
ray particles for which the correlation time can be arbitrarily 



long. In their analysis, this fact makes the time development 
of the distribution function non-Markovian , thereby precluding 
its description by a Fokker-Planck kinetic equation. They 
conclude, therefore, that the quasilinear kinetic equation 
must be dealt with in its full integro-dif f erential form. 

It is well known in plasma physics that resonance 
behavior can lead to terms in the kinetic equation that grow 
secularly with time, thus causing the quasilinear scheme to 
break down (Davidson 1972) . Specifically, the breakdown is 
due to particles whose motion is not properly described in 
the approximation scheme; over the correlation time of the 
fluctuations, orbits of such particles deviate significantly 
from the trajectories in the average field . It was our 
suspicion that quasilinear theory would fail for just this 
reason in the case considered by Klimas and Sandri (1971) . 

To test this hypothesis we have investigated a model 
equation designed to include the relevant features of the 
Klimas and Sandri model, without all of the latter's mathe- 
matical complexity. We were specifically interested in comparing 
the behavior of higher order terms in the perturbation series 
with that of the quasilinear term. We found that the quasilinear 
term in our model has the same long time behavior as that of 
Klimas and Sandri, and that the non-vanishing term of next 
higher order, while important for even fairly short times, 
strongly dominates the quasilinear term for long times. 
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We are thus led to the conclusion that in the case of 
the cosmic ray problem, the behavior of resonant particles 
cannot be adequately handled by quasilinear theory, or indeed 
by any perturbation technique which assumes that the motion of 
particles deviates only slightly from the trajectories in 
the average field. 

In §11 we describe the details of our physical model, 
emphasizing its essential similarity to that of Klimas and 
Sandri. The continuity equation in the appropriate phase 
space is used in §111 to derive a kinetic equation for the 
distribution function correct to fourth order in the strength 
of the random field, and in §IV we compare the fourth-order 
term with the second-order (quasilinear) term. An example 
with parameter values characteristic of the interplanetary 
region is discussed to show how quickly the higher order 
term dominates. 

II. The Physical Model < 

We consider a two-dimensional distribution of monoenergetic 
particles free to move in the y-z plane. A static random 
magnetic field, 6B, is oriented normal to this plane in the 
x-direction. The field is assumed to be a homogeneous 
stochastic function of z only, with Gaussian statistics and 
vanishing mean. Since any x-component of velocity as well 
as the magnitude of the velocity in the y-z plane are independent 
constants of the motion, extension of the model to three 
dimensional motion and a distribution of kinetic energies is 
trivial. Weak gradients in particle density in the z-direction 
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are allowed, but the density is assumed to be uniform in the 

y-direction. The single-particle distribution function, f, then 

can be expressed as a function b£ time , t; position, z; and 

the angle, 6, between the velocity vector, v, and the z-axis. 

fhe magnitude of v enters the analysis only as a parameter. 

We emphasize that the model, while obviously idealized, is 

nonetheless physically realizable. 

The essential feature of the Klimas and Sandri model is 
the existence of particles whose unperturbed orbits either 
remain within or periodically re-enter for arbitrarily long 
times a region throughout which values of the field strength 
are correlated. Such a class of particles, namely those for 
which 0 58 -it/2, is present also in our model. We refer to them 
as resonant because the effect of the perturbing field adds 
coherently for arbitrarily long times. 

Further similarities between the two models will be 
pointed out during the derivation of the kinetic equation for f 
in §III. 


III. The Kinetic Equation 

The continuity equation for f in the (z,0) phase space is 


9f 3 
9 1 9 z 


f > 


!_ <*£ f) 

90 l dt r; 



With dz/dt = v cos 0 and 


( 1 ) 


4| = q* B(z > g{a>( Z ), 

dt ymc ' ' 

where q is the particle charge and 
time dilation factor, Equation (1) 


( 2 ) 

y is the usual relativistic 
can be written 
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0 


(3) 


3f , « 3f ^ x 3f 

_ + VCOS0 ^ = 

We next write Equation (3) in terms of the dimensionless 

variables t = vt/z , £ = z/z and 6<f> = z 6 cjj/v, z being the 

c c c c 

correlation length of the fluctuating field: 

— + cose 35 + ae = 0 (4) 

Consider a statistical ensemble of systems, and denote 
the ensemble average of a quantity A by <A> . The fluctuating 
part of A is 6A = A - <A> ; clearly 6 <A> = 0. With this 
convention Equation (4) can be written 


3 < f ^ . Jch x u Sch 

~TT~ + cose T — + IT 

, 36f . „ ~ 36f , 3Sf 
+ IT + cos0 Tc + 6cp ~TQ 

The ensemble average of Equation (5) is 


3<f> i a 3<f> _ 
__ + COS 0 - 


<s * le > 


Subtracting Equation (6) from Equation (5) gives 


36f , A 3 6f 
IT + cos0 TT =_6(|) 


3<f> 3 6f 3 5f ^ 

30 " 30 + 30 


At this point a quasilinear analysis proceeds by neglecting 

the last two terms on the right side of Equation (7) . This 

approximation assumes that <5f is proportional to 6<f> and that 

6<j> is small compared to unity. We note that 6<p = z /r , where 

c g 

r is the gyro-radius of a particle in the fluctuating field, 

y 

Thus, the fundamental assumption of the quasilinear analysis, as 


Klimas and Sandri point out, is that z << r . Corrections to 

c g 
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the quasilinear treatment are obtained by iterating Equation 
(7) with respect to the 0(6<j> 2 ) terms. To 0(6<j> n ), 6f is thus 
written as 

5f = 6f (l) +6f (z) +... + 6f (n) (8) 

( n ) 

where f v 1 is the n-th iterate and is thus of n-th order 
in <$<(>. When the series, Equation (8), is substituted into 
Equation (6) the kinetic equation for <f>, correct through 
0 ( 5(J) n ‘ + ^) , is obtained (Kaufman, 1968). 

We carry out the iterative solution of Equation (7) via 
the method of characteristics. The characteristic unperturbed 
trajectory is defined by 


S*(t") = £*(t) + cosQ (t'-t) , 6*(x ') = 0 ,t'<x (9) 

with ?*(t) = £. The first three iterates of Equation (7), 
neglecting initial value terms, are 


6f (l) (?,0,t)^ - J dij 64) (1) 


9<f> (1) 

90 


( 10 ) 


6f U) (c,0,T)= / T dx 6<M1) / 1 dx 2 64) (2) 3 g-Q > - (2) 

o 1 o'' 0 


” f X d Ti <6^Cl) J Tl dx 6 4> ( 2 ) > (11) 

o 2 


if (s) l?,e,T >=-/ aT iS *a)^. 


X X 

+ J T dXjfiiMl)^ / 1 dx ^ <64) (2) |g. J 2 dx 3 6<M3)> 


r\ X r\ X 

+ r cix i <6(J) (1) 3 Q- / 1 dx 2 6cf> ( 2 ) / 2 dx 

o 


90- J ^ 1 3 5(j) (3 ) > 
o 


9 < 


<f> (3 ) 
90 


9<f> (3) 
90 


f> (3) 

90 


( 12 ) 
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In obtaining Equation (12) , the fact that <6<j>> = 0 was 
used. On the right side of Equations (10) - (12) an argument, 
n, means that the function is to be evaluated at the phase 
space point [C*(x n ), 0*( T n )l> a t time x^, where [?*, 9*1 is 
obtained from Equation (9) with x =x n and x = T n _i* Also, 

v x - 

Substitution of Equations (8) - (12) in Equation (6) 
gives the kinetic equation for <f> correct to fourth order 
in 6<f>. This equation can be considerably simplified by 
exploiting the statistical properties we have assumed for 6(j>, 
in particular 


a) 

< 6 cj) > = 0 , 


b) 

<6<j>(C) 54> U') 54) 

> = 


c) )> = <S(Mc)6cMr)><6<Mr') l5 <MS~ 

+ <6<j) (?) S4> U '')> <Scj) u') 6<|> u"')> + <6<H?)54>(C )> <6 4>U 

Properties (b) and (c) are true for any homogeneous Gaussian 
process for which (a) is true. Properties (a) and (b) imply 
that 6f^ 2 ^ and the third term of 6f^ 3 ^ [Equation (12) j do 
not contribute to the right side of Equation (6) . Use of (c) 
allows Equation (6) finally to be put in the form 


3<f> i — q 9<f> = JL f T 

IT 9 3? 36 [ 

+ J T dT^a+d)^ / T ' 

0 0 

+ r «t ,<«♦(! 

0 u 


di!<6* 6<t» (1)> 

dT «*(2J> la J T2 a T ,<6*6<K3J:£|f^ ) 
0 

dx 2 < 5 d> 6 4> (2)> J T2 dX 3 5d)(3)> 9 3f (3 - } ' 

o 


i )> 

) 84> (?")>■ 


(13) 
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Our contention is that in the cosmic ray problem, effects 
of <5f^ 3 ^ [the second and third terms on the right side of 
Equation (13)'] dominate those of the quasilinear term, 6f^, 
and the expansion procedure using the characteristic trajectories. 
Equation (9) , becomes invalid over times of interest. 

To facilitate evaluation of the integrals appearing in 
Equation (13) we assume that the two-point correlation 
function has the simple form 

<6<MO6<|>(c0> = < 6 2 > S(|c-e-|-l) (14) 


where S is the unit step function 


S (a-b) 



1 , a < b 
0 , a > b 


(15) 


(Recall that all lengths have been scaled by the correlation 
length, z c ) . Further following Klimas and Sandri (1971) , we 
expand <f> in a series of Legendre polynomials 


<f> = <f>o + <f> cos8 (16) 

<f> 0 is the density averaged over 8, <f> x measures the 
anisotropy and therefore is related to .the particle flux, 
and terms of higher order than P x have been neglected. Finally, 
Equations (14) - (16) are substituted in Equation (13) and 
moments of the resulting equation taken with respect to P = 1 
and P j = cos 6 to obtain coupled equations for <f> Q and < f > 1 . 

In evaluating the right side of Equation (13) it is 
assumed that spatial gradients in < f > are negligible. 
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The result is 


9< f> o + 1 9< f >i * 0 

3 x 2 3? 

1 3 <f> x 1 3<f>a 

2 3t 2 3C 


= - <6(p 2 > J dx"[I (t ";)-<6<J> 2 >I 4 (x r ) ] <f> 1 (t-t") 


where 

I 2 (t) = 


2 i y 2 

I S (T- 1 ) + i- [ (T + sin " 1 (^)]S( 1 - T ) 


(17) 


(18) 


(19) 


I 4 (x)=k^x) +k 1 (T) S(l-x) +k 2 (T) S (2-x ) (20) 

and 

V T) = ? t2 . 

kj (x) = | {( |t 2 - ~) (x 2 -l) 7 -(|x 2 -6) cos ^i) 

i 

+ 10T ^n [T+ (x 2 -l)?] } , 

k 2 <T)=i ( (r l >^ +(jT 2 -6)cos“(|) 

- 10T*n[y + (j 2 - 1) *] ), (21) 

The degree to which our simple model is successful in 

simulating the more complicated Klimas and Sandri model is 

apparent from Equations (18) and (19) . As mentioned earlier, 

our expansion parameter <Scf> 2 > = (z /r ) 2 is identical to theirs. 

c g 

Furthermore, the second-order contribution to our kernel has 
the long-time behavior I (x)^x 1 , which is exactly the long- 
time form of the quasilinear kernel obtained by Klimas and 
Sandri. Jti-st as in their model, the long range of I 2 is 
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caused by resonant particles, for which the approximation 
of unperturbed orbits breaks down. 


IV. Failure of Quasilinear Theory 

For large values of x the second and fourth-order 
contributions to the integrand on the right side of Equation 
<18) have, respectively, the forms 


<6<j> 2 > 


oo 

I 2 (t)+ I 2 (T) 



I 4 (T)-> <6(j) 2 > l“ (T) 


< 6 2 > _ 3 
12 IT T ' 


T>> 1 , 
T >>1 . 


I 4 and I 4 are plotted as functions of the dimensionless 

time variable, x, in Figure 1. From the graph it is clear 

00 

that I 4 (x) provides a lower bound for I 4 (x) for all values 
of x. It is simply verified that 1^ (x) . similarly provides 
an u PP er bound for I ^ (x) . Therefore, the ratio of the fourth 
and second-order contributions is 


I 4 (t) 


I 4 (x) 


I (x) <(S<}) > > I"(x) <<5<t) > " 24 <6( * >2> 


This ratio becomes appreciable — and the quasilinear approximation 
consequently begins to break down - on the time scale x = ( 24/<<5<J> 2 > ) 1 , 
which is comparable to or shorter than what would be the natural 
time scale of Equation (18) if I were neglected , x- <64> 2 > ~ 1 . There 
is every reason to believe that higher order terms in the 
perturbation series will become important for times < x , so that 
quasilinear theory applies, if at all, only for rather short 
times . 
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t 4 depends, of course, on the characteristics of the 
random field and the energy of the particles. For protons 
in the interplanetary magnetic field, with z^ - 2 x 10 11 cm 
and <SB 2 > 1/2 * 2.5 x 10 5 gauss (Jokipii 1971), 


< 6 cj> 2 > 


2.5 

(T+l) 2 -l 


( 22 ) 


where T is the particle kinetic energy in GeV. Since <6<j> 2 > 
must be small compared to unity for the perturbation series 
to be sensible. Equation (22) implies that the perturbation 
solution can be valid only in the energy range T > 1 GeV. For 
T = 10 GeV, one finds that = 6, corresponding to a time 
t = t z/v -40 sec. This is an order of magnitude shorter 

4 4 C 

than the time scale << 54 > 2 > 1 z c /v ^ 350 sec, illustrating 

our contention that the validity of the quasilinear approximation 

is restricted to time scales much shorter than those of interest. 
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Figure Caption 

Fig. 1. Time dependence of fourth-order contribution 

to the right side of equation (18) , both the exact 
expression, I 4 , and the approximation for long 
times, I® = t 3 /12t r. 
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